clear all
T=1;
num=50; tspan = linspace(0,T,num+1);
inits=1; options=[];
a=1; m=1; b=1;


% [time,sigma] = ode45(@myfun1,tspan,inits,options,a,b,m);
% 
% hold on 
% plot(time,sigma)
% mysol1 =@(s) m*(1-exp(-a*s))+inits*exp(-a*s);
% plot(time,mysol1(time),'r')
% s=sigmahandle;
% plot(time,s.sigma0(time,T,a,m,inits),'k')


% [time,sigma] = ode45(@myfun2,tspan,inits,options,a,b,m,T);
% 
% hold on 
% plot(time,sigma)
% mysol2 =@(s) b/a*sqrt(T)*(1-exp(-a*s));
% plot(time,mysol2(time),'r')


[time,sigma] = ode45(@myfun3,tspan,0,options,a,b,m,T);

hold on 
plot(time,sigma)
mysol3 =@(s) b*sqrt(2/T)/(a^2+(2/T)^2)*(a*(cos(2/T*s)-exp(-a*s))+2/T*sin(2/T*s));
s=sigmahandle;
plot(time,mysol3(time),'r')
plot(time,s.sigma3(time,T,a,b),'k')
